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Abstract 

Starting with a similarity function between objects, it is possible to define a distance metric on pairs of 
objects, and more generally on probability distributions over them. These distance metrics have a deep basis 
in functional analysis, measure theory and geometric measure theory, and have a rich structure that includes 
an isometric embedding into a (possibly infinite dimensional) Hilbert space. They have recently been applied 
to numerous problems in machine learning and shape analysis. 

In this paper, we provide the first algorithmic analysis of these distance metrics. Our main contributions 
are as follows: (i) We present fast approximation algorithms for computing the kernel distance between two 
point sets y and Q that runs in near-linear time in the size of T U Q (note that an explicit calculation would 
take quadratic time), (ii) We present polynomial-time algorithms for approximately minimizing the kernel 
distance under rigid transformation; they run in time 0(n-l-poly(l/e,logn)). (iii) We provide several general 
techniques for reducing complex objects to convenient sparse representations (specifically to point sets or 
sets of points sets) which approximately preserve the kernel distance. In particular, this allows us to reduce 
problems of computing the kernel distance between various types of objects such as curves, surfaces, and 
distributions to computing the kernel distance between point sets. These take advantage of the reproducing 
kernel Hilbert space and a new relation hnking binary range spaces to continuous range spaces with bounded 
fat-shattering dimension. 
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1 Introduction 



Let : M'' X M'' — > i? be a kernel function, such as a Gaussian kernel; K{p,q) describes how similar two points 
p,q are. For point sets 9, Q we can define a similarity function k^?, Q) = SpeySqeQ^^i''^)- Then the 
kernel distance is defined as 

DkCP, Q) = \/fc(a',y) + K(Q,Q)-2<y,Q). (1.1) 

By altering the kernel K, and the weighting of elements in k, the kernel distance can capture distance between 
distributions, curves, surfaces, and even more general objects. 

Motivation. The earthmover distance (EMD) takes a metric space and two probability distributions over the 
space, and computes the amount of work needed to "transport" mass from one distribution to another. It has 

become a metric of choice in computer vision, where images are represented as intensity distributions over a 
Euclidean grid of pixels. It has also been applied to shape comparison [8], where shapes are represented by 
point clouds (discrete distributions) in space. 

While the EMD is a popular way of comparing distributions over a metric space, it is also an expensive one. 
Computing the EMD requires solving an optimal transportation problem via the Hungarian algorithm, and 
while approximations exist for restricted cases like the Euclidean plane, they are still expensive, requiring either 
quadratic time in the size of the input or achieving at best a constant factor approximation when taking linear 
time [38, 2]. Further, it is hard to index structures using the EMD for performing near-neighbor, clustering 
and other data analysis operations. Indeed, there are lower bounds on our ability to embed the EMD into 
well-behaved normed spaces [4] . 

The kernel distance has thus become an effective alternative to comparing distributions on a metric space. 
In machine learning, the kernel distance has been used to build metrics on distributions [42, 23, 41, 39, 32] 
and to learn hidden Markov models [40]. In the realm of shape analysis [45, 18, 19, 17], the kernel distance 
(referred to there as the current distance) has been used to compare shapes, whether they be point sets, curves, 
or surfaces. 

All of these methods utilize key structural features of the kernel distance. When constructed using a positive 
definite^ similarity function K, the kernel distance can be interpreted through a lifting map : M'' — > to a 
reproducing kernel Hilbert space (RKHS), "K. This lifting map cj) is isometric; the kernel distance is precisely 
the distance induced by the Hilbert space Wj({{p},{q}) = \\4>{p) - 4>(.P)\\j{y Furthermore, a point set 9 has 
an isometric lifted representation $(1") = Xlpe? "^(P) ^ single vector in 'K so 0^(9, Q) = 11^(9) — ^{_Q)\\j{. 
Moreover, by choosing an appropriately scaled basis, this becomes a simple £2 distance, so all algorithmic tools 
developed for comparing points and point sets under I2 can now be applied to distributions and shapes. 

Dealing with uncertain data provides another reason to study the kernel distance. Rather than thinking 
of K{-,-) as a similarity function, we can think of it as a way of capturing spatial uncertainty; JC(p,q) is the 
likelihood that the object claimed to be at p is actually at q. For example, setting K(p,q) = exp(— ||p — 
q|p/cr)/(V^27ia) gives us a Gaussian blur function. In such settings, the kernel distance -Df (J", Q) computes the 
symmetric difference |yAQ| between shapes with uncertainty described by If. 

Our work. We present the first algorithmic analysis of the kernel distance. Our main contributions are as 

follows: 

(i) We present fast approximation algorithms for computing the kernel distance between two point sets T 
and Q that runs in near-linear time in the size ofPuQ (note that an explicit calculation would take quadratic 
time), (ii) We present polynomial-time algorithms for approximately minimizing the kernel distance under 
rigid transformation; they run in time 0(n -I- poly(l/e, logn)). (iii) We provide several general techniques for 
reducing complex objects to convenient sparse representations (specifically to point sets or sets of points sets) 
which approximately preserve the kernel distance. In particular, this allows us to reduce problems of computing 

positive definite function generalizes the idea of a positive definite matrix; see Section 2. 
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the kernel distance between various types of objects such as curves, surfaces, and distributions to computing 
the kernel distance between point sets. 

We build these results from two core technical tools. The first is a lifting map that maps objects into a finite- 
dimensional Euclidean space while approximately preserving the kernel distance. We believe that the analysis 
of lifting maps is of independent interest; indeed, these methods are popular in machine learning [41, 39, 40] 
but (in the realm of kernels) have received less attention in algorithms. Our second technical tool is an theorem 
relating e -samples of range spaces defined with kernels to standard £■ -samples of range spaces on {0, 1} -valued 
functions. This gives a simpler algorithm than prior methods in learning theory that make use of the y-fat 
shattering dimension, and yields smaller e-samples. 

2 Preliminaries 

Definitions. For the most general case of the kernel distance (that we will consider in this paper) we associate 
a unit vector and a weighting ju(p) with every p e 9. Similarly we associate a unit vector V(p) and 
weighting v(q) with every q e Q. Then we write 



where (•, •) is the Euclidean inner product. This becomes a distance Dj^, defined through (1.1). 

When y is a curve in R'^ we let U(p) be the tangent vector at p and /i(p) — 1. When 9 is a surface in M"^ we 
let be the normal vector at p and pi(p) = 1. This can be generalized to higher order surfaces through the 
machinery of fc-forms and fc-vectors [45, 35]. 

When T is an arbitrary probability measure^ in R'^, then all f/(p) are identical unit vectors and is the 
probability of p. For discrete probability measures, described by a point set, we replace the integral with a sum 
and nip) can be used as a weight k^?, Q) = ^p^y,^q^QKip,q)iJ,ip)viq). 

From distances to metrics. When K is a s)mimetric similarity function (i.e. K(p,p) — maXq^^dK{p,q), 

K{.p,q) = K{q,p), and Kip,q) decreases as p and q become "less similar") then Djf (defined through (2.1) 
and (1.1)) is a distance function, but may not be a metric. However, when K is positive definite, then this is 
sufficient for Djf to not only be a metric'^, but also for to be of negative type [16]. 

We say that a symmetric function if : E'* x R*^ ^ M is a symmetric positive definite kernel if for any nonzero 
L2 function / it satisfies J^^^ J^^^ / iq)K(p, q)fip) dpdq > 0. The proof of Dj^ being a metric follows by con- 
sidering the reproducing kernel Hilbert space associated with such aK [5]. Moreover, Dj^ can be expressed 
very compactly in this space. The lifting map : M'^ — > J{ associated with K has the "reproducing property" 
K(p,q) = (<^(p), <^(q))3^. So by linearity of the inner product, = Jp^j, <^(p)<iju(p) can be used to re- 

trieve Djf(y, Q) = — *(Q)||5^; using the induced norm || • of CK. Observe that this defines a norm 

II* Wlk = for a shape. 

Examples. If K is the "trivial" kernel, where K(p,p) — 1 and K(p, q) — Oforp ^ q, then the distance between 
any two sets (without multiplicity) T, Q is Q) - \'PAQ\, where TAQ = T U Q \ (J n Q) is the symmetric 

difference. In general for arbitrary probability measures, (CP, Q) = ||ju~'''ll2- ^fK{p,q) = (p,q), the Euclidean 
dot product, then the resulting lifting map 4> is the identity function, and the distance between two measures 
is the Euclidean distance between their means, which is a pseudometric but not a metric. 

^We avoid the use of the term 'probabiUty distribution' as this conflicts with the notion of a (Schwarz) distribution that itself plays 
an important role in the underljdng theory. 

^Technically this is not completely correct; there are a few special cases, as we will see, where it is a pseudometric [41]. 




(2.1) 
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Gaussian properties. To simplify much of the presentation of this paper we focus on the case where the 
kernel K is the Gaussian kernel; that is K{p,q) — e""^"'^" Our techniques carry over to more general 
kernels, although the specific bounds will depend on the kernel being used. We now encapsulate some useful 
properties of Gaussian kernels in the following lemmata. When approximating K{p,q), the first allows us to 
ignore pairs of points further that ■^hln^l/y) apart, the second allows us to approximate the kernel on a grid. 

Lemma 2.1 (Bounded Tails). If\\p-q\\> y/hln^l/y) then Kip,q) < y. 

Lemma 2.2 (Lipschitz). ForS^M.^ where \\5\\ < s, for points p,q ^R*^ we have \K{p,q) — K{p,q + 5)\ <s/y/h. 

Proof. The slope for xpix) = e~^^^^ is the function ipXx) = —{2/h)xe~^^^^. ^pX^) is maximized when x — 
y/h/2, which yields 'ipX\/h/2) = -\j2lhe < l/^/h. Thus \\i){x) - \i){x + s)\ < sl\fh. And since translating by 
5 changes ||p — q|| by at most e, the lemma holds. □ 



2.1 Problem Transformations 

In prior research employing the kernel distance, ad hoc discretizations are used to convert the input objects 
(whether they be distributions, point clouds, curves or surfaces) to weighted discrete point sets. This process 
introduces error in the distance computations that usually go unaccounted for. In this subsection, we provide 
algorithms and analysis for rigorously discretizing input objects with guarantees on the resulting error. These 
algorithms, as a side benefit, provide a formal error-preserving reduction from kernel distance computation on 
curves and surfaces to the corresponding computations on discrete point sets. 

After this section, we will assume that all data sets considered 7, Q are discrete point sets of size n with 
weight functions /i : T — > M and v : Q — > R. The weights need not sum to 1 (i.e. need not be probability 
measures), nor be the same for T and Q; we will set W = max(^pgj,/i(p), ^^^q v(q)) to denote the total 
measure. This impHes (since K[p,p) = 1) that k(T,T) < W^. All our algorithms will provide approximations 
of Kiy, Q) within additive error sW^. Since without loss of generality we can always normalize so that W = 1, 
our algorithms all provide an additive error of e. We also set A = (l/h)max„^gg)yQ \\u — v|| to capture the 
normalized diameter of the data. 



Reducing orientation to weights. The kernel distance between two oriented curves or surfaces can be 
reduced to a set of distance computations on appropriately weighted point sets. We illustrate this in the case of 
surfaces in R^. The same construction will also work for curves in M'^. 

For each pointp e J" we can decompose 17 (p) = iU^ip), U2ip), U^ip)) into three fixed orthogonal components 
such as the coordinate axes {e^, e^, 63}. Now 
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K(y,Q) 



if(P,q)(t/(p),V(q)) d|u(p)dv(q) 



qeQ 



qeQ 



^iV, ^SYP^^'^'^^^^'^^^ di"(p)dv(q) 



i=l 



= E 



if (P, q)(t/i(p)Vi(q)) d|L6(p)dv(q) = J] Ki^^ %\ 



1=1 



where each p e !P; has measure /ijCp) = J7[(p)|| . When the problem specifies [7 as a unit vector in R'', this 
approach reduces to d independent problems without unit vectors. 



Reducing continuous to discrete. We now present two simple techniques (gridding and sampling) to 
reduce a continuous T to a discrete point set, incurring at most sW^ error. 
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We construct a grid Gg (of size 0((A/£)'')) on a smooth shape 7, so no point"^ p e J" is further than eVh 
from a point g G G^. Let Pg be all points in CP closer to g than any other point in G^. Each point g is assigned a 
weight iu(g) = J^^p 1 d/i(p). The correctness of this technique follows by Lemma 2.2. 

Alternatively, we can sample n = 0((l/e^)(d +log(l/5)) points at random from T. If we have not yet 
reduced the orientation information to weights, we can generate d points each with weight Ui(p). This works 
with probability at least 1 — 5 by invoking a coreset technique summarized in Theorem 5.2. 

For the remainder of the paper, we assume our input dataset ? is a weighted point set in R'' of size n. 

3 Computing the Kernel Distance I: WSPDs 

The well-separated pair decomposition (WSPD) [9, 21] is a standard data structure to approximately compute 
pairwise sums of distances in near-linear time. A consequence of Lemma 2.2 is that we can upper bound the 
error of estimating K(p,q) by a nearby pair K(p,q). Putting these observations together yields (with some 
work) an approximation for the kernel distance. Since D'^('P, Q) = kC?,?) -I- k(Q, Q) — 2k(CP, Q), the problem 
reduces to computing k(J', Q) efficiently and with an error of at most (e/4)W^. 

Two sets A and B are said to be a-separated [9] if max{diam(A), diam(B)} < aminQg^ ^g^ \ \a — b\\. Let 
AiSiB — {{x,y} |xeA,yeB} denote the set of all unordered pairs of elements formed by A and B. An a-WSPD 
of a point set P is a set of pairs W = {{Ai,Bi}, . ..,{Ag,Bg}} such that 

(i) Ai,Bi CP for alii, 

(ii) AinBi = for alii, 

(iii) disjointly {Ji=iA O = P (8> P, and 

(iv) Aj and are a-separated for all i. 

For a point set P c M'' of size n, we can construct an a-WSPD of size 0{n/a'^) in time 0(nlogn+n/a'') [21, 14]. 

We can use the WSPD construction to compute D^C?, Q) as follows. We first create an a-WSPD of T U Q 
in 0(nlogn -I- n/a'^) time. Then for each pair {Ai,Bi} we also store four sets A; y = CP n A;, A; q = Q n Aj, 
Bj y = CP n Bj, and B; q = Q n Bj. Let O; € A; and e B; be arbitrary elements, and let Dj = Hoj — bi\\. By 
construction, Dj approximates the distance between any pair of elements in A; x Bj with error at most 2aDi. 

In each pair {A;,B;}, we can compute the weight of the edges from CP to Q: 

We estimate the contribution of the edges in pair (Ai,B;) to k:(CP, Q) as 

/i(a)v(b)e-^?/'^+ 2 M(b)v(a)e-^?/'^ = W^e-^?/^ 

Since D; has error at most 2aD; for each pair of points. Lemma 2.2 bounds the error as at most W;(2aDj/ Vh). 
In order to bound the total error to [s/4)W^, we bound the error for each pair by (e/4)W; since W; = 

Iip6pIiqeQi"(P)'*'(^) = Lemma 2.1, if Dj > y/hlnil/y), then e ^i^^ < y. So for any pair with 

Di > 2y/h\n{4/s), (for a < 1/2) we can ignore, because they cannot have an effect on k:(CP, Q) of more than 
(l/4)eW;, and thus cannot have error more than (l/4)£Wj. 

Since we can ignore pairs with D; > 2^/hln(^4/s), each pair will have error at most Wj(2a(2-\//iln(4/e)/ a/H) 
= Wi{4a^/]n{4/s)). We can set this equal to {s/4)Wi and solve for a < {l/4)s/y/ln{4/e). This ensures that 
each pair with Dj < 2 \/hln(4/e) will have error less than (e /4) Wj, as desired. 

■*For distributions with decajdng but infinite tails, we can truncate to ignore tails such that the integral of the ignored area is at most 
(1 — e/2)W^ and proceed with this approach using e/2 instead of e. 
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Theorem 3.1. By building and using an ((e ln(4/ s))-WSPD, we can compute a value U in time 0(n logn + 
(n/e'*)log'*/2(l/£)), such that 

|[/-d|(t,q)| < sw^. 



4 Computing the Kernel Distance II: Approximate Feature Maps 

In this section, we describe (approximate) feature representations ^>(J') = '^p^y<p{p)^J-{p) for shapes and 
distributions that reduce the kernel distance computation to an £2 distance calculation — $(Q)||j£ in 

an RKHS, JC. This mapping immediately yields algorithms for a host of analysis problems on shapes and 
distributions, by simply applying Euclidean space algorithms to the resulting feature vectors. 

The feature map allows us to translate the kernel distance (and norm) computations into operations in a 
RKHS that take time 0(np) if IK has dimension p, rather than the brute force time O(n^). Unfortunately, 'K is in 
general infinite dimensional, including the case of the Gaussian kernel. Thus, we use dimensionality reduction 
to find an approximate mapping (p ■.R'^ (where ^(7) = XIpe? "^(p)) ^^at approximates /c(y, Q): 



2 2 Kip, q)M(p)v(q) - 2 2 ( (p), (q)) 

pePqeQ psPqeQ 



The analysis in the existing literature on approximating feature space does not directly bound the dimension 
p required for a specific error bound^. We derive bounds from two known techniques: random projections [36] 
(for shift-invariant kernels, includes Gaussians) and the Fast Gauss Transform [48, 20] (for Gaussian kernel). 
We produce three different features maps, with different bounds on the number of dimensions p depending on 
logn (n is the number of points), e (the error), 5 (the probability of failure), A (the normalized diameter of 
the points), and/or d (the ambient dimension of the data before the map). 

4.1 Random Projections Feature Space 

Rahimi and Recht [36] proposed a feature mapping that essentially applies an implicit Johnson-Lindenstrauss 
projection from IK — > R'' . The approach works for any shift invariant kernel (i.e one that can be written as 
K{p,q) = fc(p — q)). For the Gaussian kernel, k{z) — e~''^" where z e W^. Let the Fourier transform of 
fc : M'' ^ M+ is g(a)) = (2;r)"''/^e""'^"^/^. A basic result in harmonic analysis [37] is that is a kernel if and 
only if g is a measure (and after scaling, is a probability distribution). Let co be drawn randomly from the 
distribution defined by g: 

kix-y)= ( gico)e^i^'^-y) dco = eM<oM,^PM)1 

where ipu>iz) = (cos((a),z)), sin((cu,z))) are the real and imaginary components of e''^"'^^ This implies that 
('0&>(^)>'0co(y)) is an unbiased estimator of k(x — y). 

We now consider a p-dimensional feature vector 0-^ : P — > where the (2i — l)th and (2i)th coordi- 
nates are described by jU(p)i/'coi(p)/(p/2) = (2M(p)cos((a>i,z))/p,2/i(p)sin(((yj,z))/p) for some coj e T = 
{co-i, cOp/2} drawn randomly from g. Next we prove a bound on p using this construction. 

Lemma 4.1. For 0^- : J' U Q ^ with p = 0^(1/ e^) login/ 5)), with probability > 1 - 5 



^Explicit matrix versions of the Johnson-Lindenstraus lemma [24] cannot be directly applied because the source space is itself 
infinite dimensional, rather than K''. 
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Proof. We make use of the following Chernoff-Hoeffding bound. Given a set [Xi, ... ,X^} of independent ran- 
dom variables, such that |Xj - £[Xj]| < A, then for M = I]|Li-^i we can bound Pr[|M - £[M]| > a] < 
2g-2a^/(nA^) ^Ve can now bound the error of using for any pair (p, q) € J" x Q as follows: 



Pr [I (</>t(p), 0T(q)) - M(p)v(q)fc(p - q)| > e/i(p)v(q)] 
= Pr [\{(l)r(pl(PM) -Er [(0t(pX '^t(<3)>] | > e^OXq)] 



< Pr 



^-/i(p)v(q) (i/'.o,(p)>^^i(q)} 



2] -M(pXq) {'ipc,t<<P)'^ojM')) 

i P 



where the last inequality follows by A < 2maXp q(2/p)ju(p)v(q) {il'co(p)>'4^<oi'i)) — 8(2/p)ju(p)v(q) since for 
each pair of coordinates ||i/)^(p)|| < 2 for all p e T (or q e Q). By the union bound, the probability that this 
holds for all pairs of points (p,q) € J" x Q is given by 

Pr [V(p,<j)ea'xQ I (<^t(p), 0T(q)) - M(p)v(q)fc(p - q)| > sui(pMq)] < {n^)2e-P''">\ 

Setting 5 > n^2e~P^^^^ and solving for p yields that for p = 0((l/e^)log(n/5)), with probability at least 
1 — 5, for all (p,q) e y X Q we have |ju(p)v(q)fc(p - q) - (0t(^)> ^tCj)) I — ^m(p)'»^(<?)- It follows that with 
probability at least 1 — 5 

peyqeQ peTqeQ peTqeQ 

Note that the analysis of Rahimi and Recht [36] is done for unweighted point sets (i.e. n{p) = 1) and actually 
goes further, in that it yields a guarantee for any pair of points taken from a manifold M having diameter A. 
They do this by building an e-net over the domain and appljdng the above tail bounds to the e-net. We can 

adapt this trick to replace the (logn) term in p by a (d log(A/e)) term, recalling A = (l/h)maXp p'gyuQ ||p— p'||. 
This leads to the same guarantees as above with a dimension of p = 0((d/£^)log(A/£5)). 



4.2 Fast Gauss Transform Feature Space 

The above approach works by constructing features in the frequency domain. In what follows, we present an 
alternative approach that operates in the spatial domain directly. We base our analysis on the Improved Fast 
Gauss Transform (IFGT) [48] , an improvement on the Fast Gauss Transform. We start with a brief review of 
the IFGT (see the original work [48] for full details). 



IFGT feature space construction. The goal of the IFGT is to approximate k{7, Q). First we rewrite k:(CP, Q) 
as the summation SqeQ^^^) where G(q) = 'v(q)2]peg)e~"^~^" nip). Next, we approximate G(q) in two 
steps. First we rewrite 

G(q) = v(q) > |U(p)e e e , 

peP 

where the quantity is a fixed vector that is usually the centroid of 7. The first two exponential terms can 
be computed for each p and q once. Second, we approximate the remaining exponential term by its Taylor 
expansion = 2i>o fr- After a series of algebraic manipulations, the following expression emerges: 

G(q) = v(q)e-T2c„(^^) 
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where C„ is given by 



The parameter a is a multiindex, and is actually a vector a = (a^, . . . , a^j) of dimension d. The expression 
z", for z e M'^, denotes the monomial z'^^z'^^ ■ ■■z'^'^, the quantity \a\ is the total degree ^ a^, and the quantity 
a! = ni(ai!). The multiindices are sorted in graded lexicographic order, which means that a comes before a' if 
|a| < |a'|, and two multiindices of the same degree are ordered lexicographically. 

The above expression for G(q) is an exact infinite sum, and is approximated by truncating the summation at 
multiindices of total degree t — 1. Note that there are at most p = C^d~^) — O^t'^) such multiindices. We 

now construct a mapping ^ : R** ^ RP. Let = J ^iw(p)e '■^ [h") ■ '^^^ 

a peP 

and S = XlqeQ G(q) is then given by 

pePqeQ a pePqeQ 



IFGT error analysis. The error incurred by truncating the sum at degree t — 1 is given by 

Err(T) = I J^J^Kip,qMp)viq) - J^J^ {^ip), Hq)) \ <J^J^Kp)v(q)^A^' = W^^A^\ 

p&Pq&Q pePqeQ psP qeQ ^" ' 

Set sW^ = Err(T). Applying Stirling's approximation, we solve for t in log(l/e) > Tlog(T/4A^). This yields 
the bounds t = O(A^) and t = 0(log(l/£)). Thus our error bound holds for t = 0(A^ + log(l/e)). Using 
p = 0(t''), we obtain the following result. 

Lemma 4.2. There exists a mapping 4> -."PuQ^RP with p = 0(A^'' + log'^Cl/e)) so 



pefqeQ pe^qeQ 



4.3 Summary of Feature Maps 

We have developed three different bounds on the dimension required for feature maps that approximate k(5', Q) 
to within eW^. 

IFGT: p = 0(A^'^ + log''(l/e)). Lemma 4.2. Advantages: deterministic, independent of n, logarithmic depen- 
dence on 1/e. Disadvantages: pol)momial dependence on A, exponential dependence on d. 

Random-points: p = 0((l/e^)log(n/5)). Lemma 4.1. Advantages: independent of A and d. Disadvantages: 
randomized, dependent on n, polynomial dependence on 

Random-domain: p = 0((d/e^) log(A/e5)). (above) Advantages: independent of n, logarithmic dependence 
on A, polynomial dependence on d. Disadvantages: randomized, dependence on A and d, polynomial 
dependence on 1/e. 

For simplicity, we (mainly) use the Random-points based result from Lemma 4.1 in what follows. If appro- 
priate in a particular application, the other bounds may be employed. 
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Feature-based computation of D^. As before, we can decompose D|(y, Q) = k(T,T) + k:(Q, Q)-2/c(J', Q) 
and use Lemma 4.1 to approximate each of fc( J", 7), k(Q, Q), and K:(y, Q) with error 

Theorem 4.1. We can compute a value U in time 0((n/e^)log(n/5)) such that \U - D|(y,Q)| < sW^, with 
probability at least 1 — 5. 

A nearest-neighbor algorithm. The feature map does more than yield efficient algorithms for the kernel 
distance. As a representation for shapes and distributions, it allows us to solve other data analysis problems on 
shape spaces using off-the-shelf methods that apply to points in Euclidean space. As a simple example of this, 
we can combine the Random-points feature map with known results on approximate nearest-neighbor search 
in Euclidean space [3] to obtain the following result. 

Lemma 4.3. Given a collection of m point sets S = {!Pi,!P2, . . . , 7"^}, and a query surface Q, we can compute 
the c-approximate nearest neighbor to Q in Q under the kernel distance in time O^pm^^'^ +o(i)) query time using 
0(pm^'^^^'^ +o(i)) space and preprocessing. 

5 Coresets for the Kernel Distance 

The kernel norm (and distance) can be approximated in near-linear time; however, this may be excessive for 
large data sets. Rather, we extract a small subset (a coreset) S from the input 9 such that the kernel distance 
between S and ? is small. By triangle inequality, § can be used as a proxy for 7. Specifically, we extend the 
notion of e-samples for range spaces to handle non-binary range spaces defined by kernels. 

Background on range spaces. Let ^(P) denote the total weight of a set of points P, or cardinality if no 
weights are assigned. Let P c R'* be a set of points and let A be a family of subsets of P. For examples of 
A, let S denote the set of all subsets defined by containment in a ball and let £ denote the set of all subsets 
defined by containment in ellipses. We say (P,A) is a range space. Let |p(A) = ^(A)/^{P). An s-sample (or 
e-approximation) of iP,A) is a subset Q c P such that 



To create a coreset for the kernel norm, we want to generalize these notions of e-samples to non-binary 
((0, l)-valued instead of {0, l}-valued) functions, specifically to kernels. For two point sets P, Q, define jc(P, Q) = 
(l/(^(P))(l/i^(Q)) XIpep SqeQ^CP'?)' ^^'^ when we have a singleton set Q = {q} and a subset P' Q P then we 
write i<p(P',q) = (l/^(P))^p(^p/K(p,q). Let = max^ ggp K'(p,q) be the maximum value a kernel can take 
on a dataset P, which can be normalized to iC"*" = 1 . We say a subset of S c P is an e-sample of (P, K) if 



The standard notion of VC-dimenion [47] (and related notion of shattering dimension) is fundamentally 
tied to the binary ({0, 1} -valued) nature of ranges, and as such, it does not directly apply to e-samples of 
{P,K). Other researchers have defined different combinatorial dimensions that can be applied to kernels [15, 
25, 1, 46]. The best result is based on y-fat shattering dimension [25], defined for a family of (0, l)-valued 
functions 3" and a ground set P. A set Y c P is ^-/at shattered by 3" if there exists a function a : Y — > [0, 1] 
such that for all subsets Z Q Y there exists some € 5" such that for every x e Z Fz(x) > a(x) + y and 
for every x e Y \ Z F^ix) < a(x) — y. Then = ^(Y) for the largest cardinality set Y c P that can be 
y-fat shattered. Bartlett et al. [7] show that a random sample of 0((l/e^)(F^log^(F^/e) -|-log(l/5)) elements 
creates an e-sample (with probability at least 1 — 5) with respect to (P, 3") for y = f^(e). Note that the y- 
fat shattering dimension of Gaussian and other symmetric kernels in M'' is d -I- 1 (by setting a(x) = .5 for 
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all x), the same as balls B in R'^, so this gives a random-sample construction for e-samples of iP,K) of size 
Oi(d/e^Xlog'{l/e) + logil/5)). 

In this paper, we improve this result in two ways by directly relating a kernel range space (P, K) to a similar 
(binary) range space (P,A). First, this improves the random-sample bound because it uses sample-complexity 
results for binary range spaces that have been heavily optimized. Second, this allows for all deterministic 
e-sample constructions (which have no probability of failure) and can have much smaller size. 

Constructions for £-samples. Vapnik and Chervonenkis [47] showed that the complexity of e-samples is 
tied to the VC -dimension of the range space. That is, given a range space iX,A) a subset Y cz X is said to 

be shattered by A if all subsets of Z c F can be realized as Z = Y n R for R G A. Then the VC-dimension of 
a range space (X,A) is the cardinality of the largest subset Y (Z X that can be shattered by A. Vapnik and 
Chervonenkis [47] showed that if the VC-dimension of a range space iX,A) is v, then a random sample Y of 
0((l/e^)(v log(l /e) + log(l /5)) points from X is an g-sample with probability at least 1 — 5. This bound was 
improved to 0((l/e^)(v -l-logl/5)) by Talagrand [44, 26]. 

Alternatively, Matousek [28] showed that e-samples of size 0((v/e^)log(v/e)) could be constructed de- 
terministically, that is there is no probability of failure. A simpler algorithm with more thorough runtime 
analysis is presented in Chazelle and Matousek [12], which runs in 0(d)^'*|X|(l/e)^^log^(l/e) time. Smaller 
e-samples exist; in particular Matousek, Welzl, and Wernisch [31] and improved by Matousek [29] show 
that e-samples exist of size 0((l/e)2-2/(^+i^), based on a discrepancy result that says there exists a labeling 
X --X ^ {-1,+1} such that raaxR^j^^^^^^^ xM < 0(|X|i/2-i/2v iQgi/2 |^|^_ alluded by Chazelle [11] 
that if an efficient construction for such a labeling existed, then an algorithm for creating e-samples of size 
0((l/e)2"2/'^''+^^log(v/e)2"^/''+^) would follow, see also Phillips [34]. Recently, Bansal [6] provided a ran- 
domized polynomial algorithm for the entropy method, which is central in proving these existential coloring 
bounds. This leads to an algorithm that runs in time 0{\X\ ■ poly(l/e)), as claimed by Charikar et al. [10, 33]. 

An alternate approach is through the VC-dimension of the dual range space (^1,^), of (primal) range space 
(X,yi), where A is the set of subsets of ranges A defined by containing the same element of X. In our context, 
for range spaces defined by balls (M'', B) and ellipses of fixed orientation (M'', £) their dual range spaces have 
VC-dimension v = d. Matousek [30] shows that a technique of matching with low-crossing number [13] along 
with Haussler's packing lemma [22] can be used to construct a low discrepancy coloring for range spaces where 
V, the VC-dimension of the dual range space is bounded. This technique can be made deterministic and runs 
in poly(|X|) time. Invoking this technique in the Chazelle and Matousek [12] framework yields an e-sample of 
size 0((l/e)2-2/^^+^^(log(l/e))2-i/^^+i^) in 0{\X\ ■ poly(l/e)) time. Specifically we attain the following result: 

Lemma 5.1. For discrete ranges spaces iX,'B) and (X,8.)forX e M.'^ of size n, we can construct an e-sample of 
size 0((l/e)2-2/Cd+i)(log(l/e))2-i/d+i) /„ 0(n • poly(l/e)) time. 

For specific range spaces, the size of e-samples can be improved beyond the VC-dimension bound. Phillips [34] 
showed for ranges "R^ consisting of axis-aligned boxes in W^, that e-samples can be created of size 0((l/e) • 
log^'^Cl/e)). This can be generalized to ranges defined by k predefined normal directions of size 0((l/e) • 
log^'^Cl/e)). These algorithms run in time 0(|X|(l/e"')polylog(l/e)). And for intervals J over M, e-samples of 
(X, J) can be created of size 0(l/e) by sorting points and retaining every e|X|th point in the sorted order [27]. 

e-Samples for kernels. The super-level set of a kernel given one input q e M'^ and a value v e M"*", is the set 
of all points p e M'' such that K{p,q) > v. We say that a kernel is linked to a range space (M.'^,A) if for every 
possible input point q € R"^ and any value v e R"*" that the super-level set of liC(-,q) defined by v is equal to 
some H G A. For instance multi-variate Gaussian kernels with no skew are linked to (R^, !B) since all super- 
level sets are balls, and multi-variate Gaussian kernels with non-trivial covariance are linked to (R'^,£) since 
all super-level sets are ellipses. 

Theorem 5.1. For any kernel iC : M x M ^ R+ linked to a range space {_M,A\ an e-sample S of {P, A) for S QM 
is a e-sample of (P, K). 
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A (flawed) attempt at a proof may proceed by considering a series of approximate level-sets, within which each 
point has about the same function value. Since S is an e-sample of {P,A), we can guarantee the density of S 
and P in each level set is off by at most 2e. However, the sum of absolute error over all approximate level-sets 
is approximately sK'^ times the number of approximate level sets. This analysis fails because it allows error to 
accumulate; however, a more careful application of the e-sample property shows it cannot. A correct analysis 
follows using a charging scheme which prevents the error from accumulating. 

Proof. We can sort all pj e P in similarity to q so that < pj (and by notation i < j) if K{pj,q) > K{pj,q). 
Thus any super-level set containing pj also contains pi for i < j. We can now consider the one-dimensional 
problem on this sorted order from q. 

We now count the deviation E{P,S,q) — kp{P,q) — i<s(S,q) from pi to p„ using a charging scheme. That is 
each element sj e S is charged to ^(P)/^{S) points in P. For simplicity we will assume that k = (^(P)/^(S) is 
an integer, otherwise we can allow fractional charges. We now construct a partition of P slightly differently, for 
positive and negative E{P,S,q) values, corresponding to undercounts and overcounts, respectively. 

Undercount of /Cs(S,q). For undercounts, we partition P into 2.^(5) (possibly empty) sets {P(,Pi,P2, 
P2,. . . ,P'^^gyP^(s)} of consecutive points by the sorted order from q. Starting with pi (the closest point to 
q) we place points in sets P- or Pj following their sorted order Recursively on j and i, starting at j = 1 and 
i = 1, we place each p; in Pj as long as K(p^,q) > K{Sj,q) (this may be empty). Then we place the next k 
points Pi into Pj. After k points are placed in Pj, we begin with Pj_^_i, until all of P has been placed in some set. 
Let t < ^(S) be the index of the last set Pj such that f (P^) = fc. Note that for all p; e Pj (for ; < t) we have 
KiSj,q) > Kipi,q), thus i<s{{sj],q) > i<piPj,q). 

We can now bound the undercount as 

?(S) ?(S) t+l 

E(P,S,q) = 2 {KpiPj,q) - i<si{Sj},q)) +Yi^p(Pj,q) <Y,^piPj,q) 
j=i j=i j=i 

since the first term is at most and since ^{Pj) — for j > t Now consider a super-level set H ^ A 
containing all points before Sf+iJ H is the smallest range that contains every non-empty P^'. Because (for j < t) 
each set Pj can be charged to Sj, then Y!i]=i n H) = fc • <^(S n H). And because S is an e-sample of (P, A), 
then ^'jtl ?(Pp = (i;;;; ?(Pp -H ?(P,- n H)) - fc • ?(S n H) < e?(P). We can now bound 

t+l t+l >, J t+l J 

Overcount of Ks(S,q): The analysis for overcounts is similar to undercounts, but we construct the partition 
in reverse and the leftover after the charging is not quite as clean to analyze. For overcounts, we partition P 
into 2^(S) (possibly empty) sets {Pj, P^, P2, P2, ■ ■ ■ , P^{s)> ^^(s)^ '^^ consecutive points by the sorted order from q. 
Starting with p„ (the furthest point from q) we place points in sets P'. or Pj following their reverse-sorted order 
Recursively on ; and i, starting at ; = ^(S) and i = n, we place each pj in Pj as long as K(pi,q) < K{Sj,q) (this 
may be empty). Then we place the next fc points pi into Pj. After fc points are placed in Pj, we begin with Pj_i, 
until all of P has been placed in some set. Let t < f (S) be the index of the last set Pj such that ?(Pp = fc (the 
smallest such j). Note that for all p; G Pj (for ; > t) we have K{sj,q) < K(j)i,q), thus i<s(.{sj},q) < i<piPj,q). 

We can now bound the (negative) undercount as 

t 1 f(S) 

£(P,S,q)= 2 {i<p(.Pj,q)-i<s({Sj},q))+ ^ ('^p(^j.q) - '^s({Sj},q)) + ^ jcpCP^q) 
7=?(S) j=t-l j=l 

1 

> (Kp(Pt_i,q)-Ks({5t_i},q)) - 2 i<s({Sj],q), 

j=t-2 
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since the first full term is at least 0, as is each Kp[Pj,q) and Kp(^Pj,q) term in the second and third terms. We 
will need the one term kp(P^_i,q) related to P in the case when 1 < (^(Pf_i) < k. 

Now, using that S is an e-sample of iP,A), we will derive a bound on t, and more importantly (t — 2). We 
consider the maximal super-level set H e >l such that no points H nP are in Pj for any j. This is the largest 
set where each point p G P can be charged to a point s e S such that Kip,q) > K[s,q), and thus presents the 
smallest (negative) undercount. In this case, HDP — U^^^Pj for somes and HnS — U^.^-^{sj}. Since t <s, then 
C(H n P) = (s - t + l)k + ^(Pt-i) = (s - t + 1)?(P)/?(S) + ?(Pt_i) and ?(H nS) = s. Thus 

^rr^^c^ fr^oD^ ' (5-t + l)g(P)/g(S) |(P^ t-l g(P,_i) 
.>?s(HnS)-?,(HnP) = — — __>_____. 

Thus (t - 2) < £?(S) + ?(Pf _i)(?(S)/?(P)) - 1. Letting = min;.ep^_^ iCCp,., q) (note K(pi, q) > q)) 

K-(Pf_i,q) K(St_i,q) ?(S) A 

EiP,S,q)> '''r - em + ^Pt-i)- 



?(P) ^(S) V J ?(S) 

, k-Kis,_^,q)-K(iP,_„q) 
— — EK -r K 



Corollary 5.1. For a Gaussian kernel, any e-sample S for (P, !B) (or for (P, £) if we consider covariance) guarantees 
that for any query q e M'^ that |fCp(P,q) - fCs(S,q)| < sK'^. 

Coreset-based computation of kernel distance. For convenience here we assume that our kernel has 
been normalized so = 1. Let P be an e-sample of CPyK), and all points p ^ P have uniform weights so 
C(P) = = W. Then for any q e M'' 

K(P,q) (c(T,q) 
s>\i<piP,q)-K^iZq)\= J^-^ ■ 

and hence 

|)c(P,q) - )c(y,q)| < = eW. 

It follows that if Q is also an e-sample for (Q,K"), then ||k:(P,Q) — kCT, Q)|| < IsW'^. Hence, via known con- 
structions of e-samples randomized [47, 44] or deterministic [28, 12, 31, 29, 43, 34] (which can be applied to 
weighted point sets to create unweighted ones [28]) and Theorem 5.1 we have the following theorems. The 
first shows how to construct a coreset with respect to D^. 

Theorem 5.2. Consider any kernel K linked with (M.'^,'B) and objects [P, Q c M c M'', each with total weight W, 

and for constant d. We can construct sets P CP and Q c Q such that |Djf (T, Q) - Dj^iP, Q)| < eW^ of size: 

• 0((l/e2-i/(d+i))iog2-i/d+i(i/g))^ via Lemma 5.1; or 

• 0(l/e^)(d -I- log(l/5))) via random sampling (correct with probability at least (1 — 5)). 

We present an alternative sampling condition to Theorem 5.2 in Appendix A. It has larger dependence on e, 
and also has either dependence on A or on log n (and is independent of iC"'"). Also in Appendix A we show that 
it is NP-hard to optimize s with a fixed subset size k. 

The associated runtimes are captured in the following algorithmic theorem. 

Theorem 5.3. When K is linked to (R'^, !B), we can compute a number D such that |I>jf (CP, Q) - b\ < s in time: 

• 0(n • (l/e2d+2)iogd+i(i/g)); or 

• 0(n -I- (logn) • ((l/e^)log(l/5)) -I- (l/e'^)log^(l/5)) that is correct with probability at least 1-5. 

Notice that these results automatically work for any kernel linked with {M.'^,'B) (or more generally with 
(M'', £)) with no extra work; this includes not only Gaussians (with non-trivial covariance), but any other 
standard kernel such as triangle, ball or Epanechnikov. 
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6 Minimizing the Kernel Distance under Translation and Rotation 



We attack the problem of minimizing the kernel distance between CP and Q under a set of transformations: 
translations or translations and rotations. A translation T e M'* is a vector so Q ® T = {q + T | q e Q}. The 
translation T* = argmin^gj^d Djf(T, Q © T), applied to Q minimizes the kernel norm. A rotation R e SO(d) can 
be represented as a special orthogonal matrix. We can write RoQ = {i?(q) | q e Q}, where i?(q) rotates q about 
the origin, preserving its norm. The set of a translation and rotation (T*,R*) = argmin(^j-jj)g]gdxSO(d) ■^^icC?''^ ° 
(Q ® r)) applied to Q minimizes the kernel norm. 

Decomposition. In minimizing Djf(CP,i?o (Q © T)) under all translations and rotations, we can reduce this 
to a simpler problem. The first term fc([P, J") = Xlpjsy Sp263''"^Pi-^'^'^^^2)^(Pi>P2) has no dependence on T or 
R, so it can be ignored. And the second term k(Q, Q) = '^q^^Q'^q^^Qv(^qi)v{q2)K{R{qi + T),R(q2 + T)) can 
also be ignored because it is invariant under the choice of T and R. Each subterm K"(i?(qi + T),R(q2 + T)) only 
depends on \ \R(,qi + T) — R(q2 + T)!! = Hq^ — q2\\, which is also independent of T and R. Thus we can rephrase 
the objective as finding 

(T*,r) = arg max V yM(p)v(q)K(p,i?(q + T)) = arg max k('?,Ro(Q®T)). 
(r,R)eR'ixSO(d)^^ (r,R)e]RdxS0(d) 

We start by providing an approximately optimal translation. Then we adapt this algorithm to handle both 
translations and rotations. 

6.1 Approximately Optimal Translations 

We describe, for any parameter e > 0, an algorithm for a translation t such that D^(CP, Q © T) — d|^(CP, Q © T* ) < 
sW-^. We begin with a key lemma providing analj^ic structure to our problem. 

Lemma 6.1. )c(y, Q © T*) > W^/n^. 

Proof. When T e M'' aligns q G Q so q + T = p for p G P it ensures that K(p,q) = 1. We can choose 
the points p and q such that ju(p) and v(q) are as large as possible. They must each be at least W/n, so 
K(p, q)/i(p)v(q) > W^/n^. All other subterms in k(T, Q © T) are at least 0. Thus k(T, Q © P) > W^/n^. □ 

Thus, if k(T, Q © T*) > W^/n^, then for some pair of points p e CP and q e Q we must have jLi(p)v(q)iC(p, q + 
T*) > ju(p)v(q)/n2, i.e. iC(p,q + P*) > l/n^. Otherwise, if all n^ pairs (p,q) satisfy ju(p)v(q)iC(p,q + P*) < 
M(pXq)/n^, then 

k(,Z Q © P*) = Y,Y,l^(PMqmp,q + P*) < XiZiM(p)v(q)/n2 ^ p^2/„2_ 

Thus some pair p e CP and q e Q must satisfy ||p — (q + P*)|| < -\/ln(n^, via Lemma 2.1 with y — l/("^)- 

Let Gg be a grid on R'^ so that when any point p e M'^ is snapped to the nearest grid point g e Gg, we 
guarantee that ||g — p|| < s. We can define an orthogonal grid Gg — {{slyfd)z \ z e Z'^}, where is the 
d-dimensional lattice of integers. Let S[£,p,A] represent the subset of the grid Gg that is within a distance A 
of the point p. In other words, g[e,p,A] = {g € Gg | ||g -p|| < A}. 

Algorithm. These results imply the following algorithm. For each point p e CP, for each q e Q, and for each 
g e 5[e/2,p, \/ln(n^)] we consider the translation Tp q ^ such that q + Tp ^ ^ = g. We return the translation 
Tp q g which maximizes fc(CP, Q © Tp q g), by evaluating k at each such translation of Q. 

Theorem 6.1. The above algorithm runs in time 0{(1/ eYn'^ log'^^^ n), for fixed d, and is guaranteed to find a 
ti-anslation t such that dU?, Q © P) - D|(CP, Q © P*) < eW^. 
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Proof. We know that the optimal translation T* must result in some pair of points p e ^ and q e Q such 

that Hp — (q + T*)|| < ^/ln(n^) by Lemma 6.1. So checking all pairs p G T and q G Q, one must have 

I |p - q 1 1 < \/ln(n2). Assuming we have found this closest pair, p e CP and q e Q, we only need to search in the 
neighborhood of translations T = p — q. 

Furthermore, for some translation Tp ^ ^ = g - q we can claim that fc(CP, Q © T*) - k:( CP, Q © g g) < e. Since 
\\T* - Tp q gW < s/2, we have the bound on subterm |A:(p,q + T*)-K{p,q + ^ g)| < e/2, by Lemma 2.2. In 
fact, for every other pair p' e T and q' e Q, we also know \Kip',q' + T*) - K[p',q' + Tp q g)\ < s/2. Thus the 
sum of these subterms has error at most (£/2)^pg,j,^ggg/i(p)v(q) = (e/2)W^. 

Since, the first two terms of D^('P,Q®T) are unaffected by the choice of T, this provides an e-approximation 
for D^C?, Q © T) because all error is in the (— 2)fc(CP, Q © T) term. 

For the runtime we need to bound the number of pairs from CP and Q (i.e. 0(n^)), the time to calculate 
k:(CP, Q© r) (i.e. 0(n^)), and finally the number of grid points in 9[e/2,p, V^lnC"^]- The last term re- 
quires 0((l/sY) points per unit volume, and a ball of radius ^/ln(rL^) has volume 0(log''^^n), resulting in 
O^il/sY log'^/^ n) points. This product produces a total runtime of Oiil/sYn'^ log'^^^ n). □ 

For a constant dimension d, using Theorem 5.2 to construct a coreset, we can first set n = 0((1/ e^) log(l /5)) 
and now the time to calculate k^'J', Q © T) is 0((l/e'^)log^(l/5)) after spending 0(n + (l/e^)log(l/5)logn) 
time to construct the coresets. Hence the total runtime is 

0(n + logn(l/e2)(log(l/5)) + (l/e-^+S) • log'^/2((l/8)log(l/5))log4(l/5)), 
and is correct with probability at least 1 — 5. 
Theorem 6.2. For fixed d, in 

0(n + logn(l/e2)(log(l/5)) + (1/8^+^) ■ log'^/'((l/s) log(l/5))log4(l/5)) 
time we can find a translation f such that 0^(7, Q © f ) - DfCCP, Q © T*) < sW^, with probability at least 1-5. 

6.2 Approximately Optimal Translations and Rotations 

For any parameter e > 0, we describe an algorithm to find a translation t and a rotation R such that 

(CP,R o (Q © f )) - dI(J>,R* o (Q © T*)) < eW^. 

We first find a translation to align a pair of points p e J" and q e Q within some tolerance (using a method 
similar to above, and using Lemma 2.1 to ignore far-away pairs) and then rotate Q around q. This deviates 
from our restriction above that R e SO(d) rotates about the origin, but can be easily overcome by performing 
the same rotation about the origin, and then translating Q again so q is at the desired location. Thus, after 
choosing a q e Q (we will in turn choose each q' e Q) we let all rotations be about q and ignore the extra 
modifications needed to f and R to ensure R is about the origin. 

Given a subset S c Q of fewer than d points and a pair (p, q) e CP x Q where q ^ S, we can define a rotational 
grid around p, with respect to q, so that S is fixed. Let 31^ g be the subset of rotations in d-space under which the 
set S is invariant. That is for any R^Jl^s ^^d any s e S we have R{s) = s. LetT = d — |S|. Then (topologically) 
^d,s = SO(t). Let Rs^p^q = minKe3?^j ||i?(q) -p|| and let q = Rs,p,q(.q)- Let Ji[p,q,S,s,A] c Di^ s be a set 
of rotations under which S is invariant with the following property: for any point q' such that there exists a 
rotation R' e Jl^s where R'iq) = q' and where \\q' — q\\ < A, then there exists a rotation R G 'K{p,q,S,e,K\ 
such that \ \R{q) — q!\\< e. For the sanity of the reader, we will not give a technical construction, but just note 
that it is possible to construct CK[p,q,S,e, A] of size 0((A/e)'^). 
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Algorithm. For each pair of ordered sets of d points (pi,P2, • • • ,Pti) c ? and (qi, q2. • • • . ^d) c Q consider the 

following set of translations and rotations. Points in T may be repeated. For each g e S[e/d,p^, -yin(max{l/£', n^})] 
consider translation Tp^ ^^ g such that + Tp^ q_^ g = g. We now consider rotations of the set (Q © Tp_^ q^ g). 

Let S = {qi} and consider the rotational grid Ji[p2,q2 + '^Pi,qi,S'^'^l^' VlnCl/e)]. For each rotation R2 ^ 

^>^[P2,q2 + "^Pi.qi.g'S.e/d, i/ln(l/e)] we recurse as follows. Apply i?2(Q © ^pi,qi,g) and place ^2(52 + ^pi.gi.g) 

in S. Then in the ith stage consider the rotational grid 'K[pi,Ri_i{. . .i^2(q2 + ^pi,qi,g) • • O^S, e/d, -^/inCl/g)]. 
Where is some rotation we consider from the ith level rotational grid, let R — R^o R^_i o . . . o R^. Let {f,R) 
be the pair (Tp g that maximize /c(?,Ro (Q © Tp g g)). 

Theorem 6.3. The above algorithm runs in time 

0(n2''+2^/^)(d^-d+2)/2iQg(d^-3d+2)/4(i/^)l0gd/2(j^^^^„2^;^/^j))^ 

for a fixed d, and is guaranteed to find a translation and rotation pair (^f,R), such that 

Dli9,Ro (Q © f )) - Dli9,R* o Q © T*) < sW^. 



Proof. We compare our solution {_f,R) to the optimal solution (T*,i?*). Note that only pairs of points (p,q) € 
y X Q such that \ \p-R*(q + T*)\\ < ^/ln{l/e) need to be considered. 

We first assume that for the ordered sets of d points we consider (pi,P2j--- jPd) ^ ^ i<ii><l2> - ■ ■ ><ld) ^ 2 
we have (Al) Hp; - R*(qf + T*)|| < v^ln(l/e), and (A2) for S - {qi, . . . ,qi_i}, let q^ e Q be the furthest point 
from S such that ||pi - (q^ + T*)|| < ^ln(l/£-). Note that (A2) implies that for any rotation R e 31^ s that 
llq; - J?(qi)|| > llq' -Riq')\ \ for all q' e Q that can be within the distance threshold under (T*,R*). In the case 
that fewer than d pairs of points are within our threshold distance, then as long as these are the first pairs 
in the ordered sequence, the algorithm works the same up to that level of recursion, and the rest does not 
matter. Finally, by Lemma 6.1 we can argue that at least one pair must be within the distance threshold for our 
transition grid. 

For each point q e Q we can show there exists some pair (r,i?) considered by the algorithm such that 
\\R*iq + T*) — R(q + T)|| < s. First, there must be some translation T — g in our grid that is within a 
distance of e/d of T*. This follows from Lemma 2.2 and similar arguments to the proof for translations. 

For each we can now show that for some i?j e 3-C (the rotational grid) we have \\R^(R^_l(. . .R2(qi + 
Tp^q^ g)...)) — R*(qj + T*)|| < e. By our assumptions, the transformed q^ must lie within the extents of "K. 
Furdiermore, there is a rotation R'j that can replace each Rj for ; e [2, i] that moves q; by at most e/d such 
that ..R'2(.qi) ...)) = R*(.qi)- Hence, the composition of these rotations affects q^ by at most — 1), 

and the sum effect of rotation and translation errors is at most e. 

Since each qi is invariant to each subsequent rotation in the recursion, we have shown that there is a pair 
{.T,R) considered so \\R{.qi + T) - R*{_qi + T*)\\ < s for q^ in the ordered set (qi,q2, • • • ,qd)- We can now use 
our second assumption (A2) that shows that at each stage of the recursion q; is the point affected most by the 
rotation. This indicates that we can use the above bound for all points q G Q, not just those in our ordered set. 

Finally, we can use Lemma 2.2 to complete the proof of correctness. Since if each K{p, q) has error at most 
£, then 



XiZi/^(i'>(^)^(P'^Cq + f)) - XiZiM(p)v(qK(p,i?*(q + T*)) 



We can bound the runtime as follows. We consider all d! (|^) = 0(n'') ordered sets of points in Q and all n'^ 
ordered sets of points from y. This gives the leading 0(n^'^) term. We then investigate all combinations of grid 
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points from each grid in the recursion. The translation grid has size 

The size of the ith rotational grid is 0{{'sJ\o%{\ / e) I eY~^ , starting at i = 2. The product of all the rota- 
tional grids is the base to the sum of their powers Xfj/C'^ ~ = Sfj/ i = {d — l)(d — 2)/2 = (d^ - 3d + 
2)/2, that is 0((l/e)'^''''"^''+2^/^ log'^'^'"^'^+^^/'*(l/e)). Multiplying by the size of the translational grid we get 
0((l/e)':'''-'*+2)/2iog(d2-3d+2)/4Q/^-)lQgd/2(-jj^^^^y^2^;^/g|)-)_ jl^gj^ g^^l^ rotation and translation we must 

evaluate K:(y,i? o (Q ® r)) in O(n^) time. Multiplying these three components gives the final bound of 

The runtime can again be reduced by first computing a coreset of size 0((l/e^)log(l/5)) and using this 
value as n. After simplifying some logarithmic terms we reach the following result. 

Theorem 6.4. For fixed d, in 

0(n + logn(l/e2)(log(l/5)) + il/sf'^'+^'^+^^'^(logil/e5)f'^'+'^'^+^°'^^'^), 
time we can find a ti-anslation and rotation pair if,R), such that 

0^(9, Ro (Q e t)) - D^(?,R* o Q e T*) < eW^, 

with probability at least 1 — 5. 
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A Coresets for the Kernel Distance 
A.l Alternative Coreset Proof 

In Section 5 we presented a construction for a coreset for the kernel distance, that depended only on 1/e for a 
fixed kernel. This assumed that K'^ — maxp ^ Kip, q) was bounded and constant. Here we present an alternative 
proof (also by random sampling) which can be made independent oiK'^. This proof is also interesting because 
it relies on p the dimension of the feature map. 
Again, we sample k points uniformly from P and reweight -qip) = W /k. 

Lemma A.I. By constructing S with size k = 0{il/ e^)login/ 5)\ogi(l/ e5)logn)) we guarantee dI-['?,§) < sW'^, 
with probability at least 1 — 5. 

Proof. The error in our approximation will come from two places: (1) the size of the sample, and (2) the 
dimension p of the feature space we perform the analysis in. 

Let 4> : y X R"*" "K describe the true feature map from a point p e y, with weight ju(p), to an infinite 
dimensional feature space. As before, set = Xlpey At(p)0(p)' recall that (J, Q) = — 
for any pair of shapes T and Q. 

By the results in the previous section, we can construct ^ : ^ x R"*" W (such as 0-f defined for Lemma 4.1) 
such that !>(?') = Spey for ^^ly pair of shapes J" and S with weights W — 2]pga)M(p) — Spes 'HCp)? 

we have |||*(a') - *(S)||^ - - l>(S)|P| < (e/2)w2, with probability at least 1 - 5/2. This bounds the 

error in the approximation of the feature space. 

We now use the low dimension p of this approximate feature space to bound the sampling error. Specifically, 
we just need to bound the probability that ||*(a') - *(S)f = ||E[*(S)] - *(S)||^ > (e/2)w2, since £[*(§)] = 
!>(?). This is always true if for each dimension (or pair of dimensions if we alter bounds by a factor 2) 

m e [l,p] we have |l'(§)m — E[l'(S)m] | < \/ eW^/(2p), so we can reduce to a 1-dimensional problem. 

We can now invoke the following Chernoff-Hoeffding bound. Given a set {X^, . . . ,Xj.} of independent random 
variables, such that -£[Xj]| < A, then for M = X;i^=i we can bound Pr[|M-r£[Xi]| > a] < 2e~'^"-''l^''^'''> . 
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By letting a — yfw^^ej{2p) and — 4>(Pi)m> the mth coordinate of ^(p;) for p; e S, 



Pr 



l>(S) - > (e/2)w2 



= Pr 

<pPr 

<p2e 



*(§)-£[!'(§)] >(e/2)w2 



where the last inequality follows because A = max^ — < 2W/k since for any p e S we have 

||0(p)|| = w/fc. By setting 5/2 > p2e~*^^/^4p)^ can solve for fc = 0((p/e)log(p/5)). The final bound 
follows using p = 0((l/e^)log(n/5)) in Lemma 4.1. □ 

Again using the feature map summarized by Lemma 4.1 we can compute the norm in feature space in O(pfc) 
time, after sampling k — 0((l/e'^)log(n/5)log((l/e5)logn)) points from ? and with p = 0((l/£^)log(n/5)). 

Theorem A.I. We can compute a vector S - ^CS) e Rf in time 0(n + (l/e^)log^(n/5)log^((l/e5)logn))) such 
that - S||^ < eW^ with probability at least 1-5. 

A. 2 NP-hardness of Optimal Coreset Construction 

In general, the problem of finding a fixed-size subset that closely approximates the kernel norm is NP-hard. 

Definition A.1 (Kernel Norm). Given a set of points 9 — {Xi}^^^, a kernel function K, parameter k and a 
threshold value t, determine if there is a subset of points S c ? such that |S| = k and Dj^ (§, IP) < t. 

Theorem A.2. Kernel Norm is NP-hard, even in the case where k — n/2 and t — 0. 

Proof. To prove this, we apply a reduction from Partition: given a set Q = {xj"^^ of integers with sum to 
2"^^ = 2m, determine if there is a subset adding to exactly m. Our reduction transforms Q into a set of points 
'9 = {^i}"=i which has subset S of size k = n/2 such that ||S - T|| < t if and only if Q has a partition of two 
subsets Qi and Q2 of size n/2 such that the sum of integers in each is m. 

Let c = j^^Xi = 2m/n and x- = Xi — c and let t = 0. Let the kernel function K be an identity kernel defined 
K(a, b) — (a, b), where the feature map is defined <^(a) = a. This defines the reduction. 

Let s = (n/fc) Xx'es P ~ Sx'eO' ^i'- Note that p = by definition. Since we have an identity kernel so 

0(a) = a, (S, J") = ||s — p||. Thus there exists an S that satisfies Dj^ (§, T) < if and only if s = 0. 

We now need to show that 5 can equal if and only if there exists a subset c Q of size n/2 such that its 
sum is m. We can write 



Thus s = if and only if 2]x'eS ~ Since S must map to a subset Qi c Q, where x^ e S implies X; e Qi, 
then s = holds if and only if there is a subset Qi c Q such that eQi ~ would define a valid 

partition, and it completes the proof □ 
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